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Abstract 

>' 
iy-\ ' Using measured instrumental response functions for data deconvolution is a known source of 

._!. , uncertainty. This problem is revisited here with Bayesian data analysis an Monte Carlo simula- 

vQ ■ tions. Noise correlation induced by the convolution operator is identified as a major source of 

^^ ' uncertainty which has been neglected in previous treatments of this problem. Application to a lu- 

O , minescence Ufetime measurement setup shows that existing approximate treatments are markedly 

^.^' defficient and that the correlation length of the noise is directly related to the lifetime to be esti- 

Oj, mated. Simple counteractive treatments are proposed to increase the accuracy of this procedure. 

> 

X 



1 Introduction 

The deconvolution problem is a classical inverse problem, and has received a lot of attention in many 
scientific and engineering fields. The instrument response function (IRF), also called blurring func- 
tion, is generally assumed to be accurately determined. It is however not uncommon that the IRF is 
measured with the same accuracy as the signal to be treated, due to instrumental or experimental de- 
sign constraints. The impact of an uncertain IRF on the accuracy of the deconvolved signal has to be 
considered with care. This uncertainty propagation issue has been addressed in the past by Dose et al. 
IjJ]. We show here analytically and numerically that their approximate solution does not encompass 
important effects of noise correlation due to convolution. 

In this paper, we use Bayesian data analysis to derive an exact expression of the likelihood function 
in the case of gaussian additive noise. This solution is applied to the classical problem of lifetime 
estimation from luminescence data. 

2 Theory 

The observed signal vector s (length n) is generally expressed as a linear reconvolution model 



s = Hm + e^ 



(1) 



where m is a vector of values of the model function at the measurement points, H is a n x n zero- 
padded lower triangular Toeplitz matrix built from the IRF h of length Uh 
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(2) 



and e^ is an additive noise with multinomial statistics and covariance matrix R^ : 

e, ~AA„(0,R,). (3) 

Note that we use bold lowercase symbols for vectors (s) and bold capitals for matrices (H). 
Using the symmetry property of convolution, Eq. ^can also be written 

s = Mh + Bs (4) 

where M is a n x n/i lower triangular Toeplitz matrix built from the model vector m as shown above 
(Eq.0. 

As the exact IRF is generally not known, a measured IRF is used instead to solve EqQ or Eq@] A 
multinormal additive noise model is also used for the IRF 

h = h + eh, (5) 

where h is the unknown exact IRF and e^ ~ J\fn,^ (0, R/i). 

The problem is to reconstruct the model vector m, knowing s, R^, h and R/j, and to evaluate the 
impact of the measurement uncertainties of h on m. This is an inverse problem doubled with an 
uncertainty propagation problem. Bayesian data analysis is very well suited to handle this kind of 
problem libl lib I- 

2.1 Bayesian data analysis 

The posterior probability density function (pdf) for m is obtained by Bayes's formula 

pirn) 
p(m|s,Rs,h,R/,) = ^^— p(s|m,Rs,h,R/,), (6) 

p(s) 

where p{m) and p{s) are the prior pdf's for m and s, and where p(s|in, Rs, h, R/i) is the likelihood 
function. Given our model, we do not know explicitely this latter function. Instead, we know the 



explicit expression for the likelihood when the exact IRF h is considered (cf . Eq. |3ll 



p(s|m,R„h) ~AA„(Mh,R, 



(7) 



Applying the marginalization rule and knowing the expression of the pdf for h 



p(h|h,RO~AA„,(h,R^), 



(8) 



we can write 



p(m|s,Rs,h,R/,) 



p(m) 
p{s) 



dh p{s\m, Rs, h)p(h|h, R/^). 



(9) 



Considering that in our model p(s) is a normalization constant, and expliciting the pdf 's, one gets 



p(m|s,Rs,h,R/i) cxp(m) / dh exp ( --J 



(10) 



where 



J = (s - Mh)^ R7Xs - Mh) 
+ (h-hfR,i(h-h). 



(11) 



This quantity is rearranged in order to enable analytical integration 



\Tt3-1 



T-D-U 



J = (h-ho)^P-^(h-ho)-h^p-^ho 



+ s^R7^s + h^R^^h, 



(12) 



where 



P = (M^R^^M + R-^)"^ 
ho = P(M^R7is + R^^h) 
= h + PM'^R7i(s-Mh) 



(13) 



Integration over h finally leads to 

p(m|s, R„ h, R;,) (X ^^ exp f-^(s - Mh)^K(s - Mh)) 



|p|i/2 "^"^ \^ 2^" """^ "'^" ""V / > (14) 

where 

K = R;1 - R^^MPM^R^^ (15) 

This expression for the posterior pdf calls for a few comments : 

• convolution of the model vector with a noisy IRF leads to a "noisy model" m = Mh, affected 
by correlated noise with covariance matrix K^^, the structure of which depends explicitely on 
the model vector itself (heteroscedastic correlated noise). This is in contrast with the result of 
Dose et al. ji]/, who obtain an expression for an effective variance, and do not consider the 
covariance part. 

• the mode of the posterior pdf depends on the actual value of the measured IRF h. A bias in the 
optimal values for the model vector is thus to be expected, as a different realization of the IRF 
would lead to a different solution. In any case, a consistent uncertainty analysis should ensure 
that the exact value lies within confidence intervals. 

3 Application 

An application of interest is for instance the lifetime estimation of unstable chemical species from 
their luminescence decays. 

3.1 Model 

A mono-exponential decay signal with lifetime r is generated over a regular time grid (n/j = n = 
100). The model is rrii = exp(— ij/r). The IRF is a gaussian function centered at t^, and of FWHM 

hi = exp (-41n(2)(ti - tof/wl) . (16) 



In order to keep a single parameters, the model after convolution is rescaled to the maximal value of the 
signal, and we can set p(t[s, R^, h, R/j) = p(m|s, R^, h, Rh). Homoscedastic noise is considered 
for both signal and IRF, i.e. R^ = a^ * I„ , R/, = cr^ * I„^. Finally, a uniform prior distribution for r 
is used (j){t) = cte). 



3 0.5 



0.2 



0.4 







- 


i \ 


signal 

-- IRF 






■^-^V^^U^ 


'y^ 



0.8 



time (arb. units) 



Figure 1: Typical synthetic signal and IRF used for lifetime estimation (o-^ 
to = O.land Wh = 0.03). 



ah = 0.01, r = 0.1 



3.2 Comparison of models of the posterior pdf 



The exact expression for the posterior pdf (eallSt is compared to approximate expressions : 



no correction for the noisy IRF (cj/j = 0, in our model), which is the most commonly used 
method; 

the diagonal approximation of eq^J which implements some level of variance correction, but 
fails to encompass the correlation in the model's noise; 



the "effective variance" method 



uiyi 



Variance. Fig. IJlrepresents the posterior pdf p(r|s, cis, h, ah) computed by Monte Carlo simulation, 
and by the different methods in the case of a same measurement accuracy for the signal and the 
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Figure 2: Posterior density functions for the lifetime estimated from a noisy decay and for different 
treatments of the IRF's uncertainty (a^ = ag = 0.01). The functions have been shifted and renor- 
malized to facilitate direct comparison, (a) comparison of simulation results (1000 runs) with the 
full bayesian solution proposed in the present work; (b) comparison of the full treatment with vari- 
ous approximations (see text), full support of the IRF; (c) support of the IRF limited to t < 0.2 (all 
approximate methods are undiscemable). 



IRF (ah = as = 0.01). The Monte Carlo method consists in repeated analysis of randomly noised 
signal and IRF to build histograms of the maximum a posteriori (MAP) lifetime values (modes of the 
posterior pdf). All curves have been shifted to a common mode, in order to facilitate comparison. The 
exact expression is fully coherent with the histogram resulting of the simulation, i.e. it takes properly 
the variance of the signal and the variance of the IRF into account. It can be seen on this figure that the 
approximate methods all perform quite similarly and fail to recover the full variance of the lifetime. 
The "effective variance" method is seen to be numerically equivalent to the diagonal approximation 
of our method, and it performs only slightly better than the totally uncorrected method. Correlation in 
the noise of the convolved model can thus have a major impact on uncertainty quantification. 

Bias. All methods perform similarly with regard to the bias on lifetime estimation (Fig. |3ll. In this 
figure, we reported the estimation by the "effective variance" method as function of the estimation by 
the exact method. The biases of both methods are highly correlated and practically identical. How- 
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Figure 3: Error estimates and 95% confidence intervals for the lifetime recovered simultaneously by 
the exact method Argxact and by the effective variance method Arejj for 100 randomly noised signals 
and IRF's {ay, = (Ts = 0.01). 

ever, underestimation of the confidence intervals by the approximate method results in inconsistent 
estimations, i.e. it fails significantly more than the exact method to include the exact value inside 
the confidence interval, and confidence intervals for different realizations of the noise are frequently 
disjoint. In this regard, Eq. [TSlperforms much better. 

Accuracy of the IRE For a given lifetime, when the IRE is measured with a better accuracy (uh < 
as), the differences observed between the various methods tend to vanish (EiglU. Eor instance, if the 
IRE is ten times more accurate than the signal, the uncorrected method provides exact results over all 
the practical range of lifetimes. It is also observed that longer lifetimes are relatively more affected 
than shorter ones, which is a pure effect of noise correlation (see next section). 

3.3 Structure of the correlation matrix 



The convolution of the mono-exponential model by the IRE is a vector lii which elements obbey the 

following reccurence 

At 

(17) 



rhj = exp( )mi-i + hj + Sh- 
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Figure 4: Ratio of the standard deviations of the posterior pdf for the exact model {a exact) and for the 
uncorrected model (sJuncm) as a function of the theoretical Ufetime. 

As soon as the IRF vanishes the correlation between consecutive points is 



< mj,mj_i >= exp( 



At, 



(18) 



As there is supposedly no correlation in the signal noise, the covariance matrix K preserves this 
correlation scheme. An approximation of the correlation matrix can thus been expressed as 
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(19) 



where p = exp(— -^). The noise correlation decays thus exponentially with de delay between points 
of the model, the decay rate being the inverse of the theoretical lifetime. 



3.4 Support length of the IRF 

When the IRF is recorded on the same support as the signal, most of the its elements are pure noise. 
We saw above that these points contribute significantly to the correlation of the noise in the convolved 
model. Limiting the support of the IRF (rih < n), or zeroing it's purely noisy elements might thus 
enable to improve the correlation matrix. If we observe the standard deviation for the convolved 
model (Fig. [Sj, we see that the truncation of the support of the IRF contributes significantly to reduce 
the uncertainty at larger times. As shown on Fig|2c), this enables some uncertainty reduction for the 
hfetime estimation, but the effect of the correlated noise is still quite marked. 



> 
-a 
-a 



0.025 


- 


1 




A-\rs 




" Aa' 


1 


.^y^\yjWy^ 




0.02 


~1 


1 
1 




~ 


0.015 




\ 

1 
1 


\ 




- 




\ 
\ 


FulllRF 

- - Truncated IRF 


0.01 


\ 

\ 




0.005 


- 


1 


\ 
\ 

\ 
\ 
\ 
\ 
\ 
\ 

1 V~-i — , 1 


- 



time (arb. units) 



Figure 5: Standard deviation of the model, after convolution with the IRF. Full line : IRF with full 
support; dashed line : IRF with support limited to t < 0.2. Same parameters as for Fig. |2l 



3.4.1 Identifiability 

The behaviour of the present model with regard to the limits of detection of lifetimes due to the IRF 
has been tested by reconstructing the posterior pdf from synthetic signals generated with very small 
lifetimes. The posterior pdf displays explicitely the non identifiability of lifetimes that are too small 
(Fig. |6ll. When r decreases, the pdf becomes asymmetric, defining an upper limit for the lifetime, but 
no lower limit, except the one imposed by the prior. 
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Figure 6: Evolution of the (unnormalized) posterior pdf for r at the resolution limit of the experimental 
setup. The exact values for r are reported alongside the curves. The standard deviation for signal and 
IRF is (js = ah = 0.005. 

4 Conclusion 



The use of measured instrumental response functions for data deconvolution is a source of uncertainty. 
We derived a new expression of the likelihood within a bayesian framework to explicitely incorporate 
this effect and display it's importance. Convolution of a noisy IRF with a model curve produces a 
noisy model curve with correlated noise. 

This has been illustrated on a luminescence lifetime measurement setup, for which it was shown that 
existing approximate treatments were markedly defficient. It was also shown that, in this case, the 
correlation length of the noise was directly related to the lifetime to be estimated. Longer lifetimes 
are thus counterintuitively more affected by IRF's uncertainty that shorter ones. Although the most 
efficient way to reduce this effect is clearly to improve the IRF's measurement accuracy, we have 
shown that an qualitative improvement can very simply be obtained by zeroing those parts of the IRF 
consisting of pure noise. 

The method has been applied to an homoscedastic noise pattern, but extension to cases where the 
noise is dependent on signal intensity (e.g. photon counting methods) is straightforward, as long 
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as the normal noise distribution approximation is valid. Similarly, cases where the IRF is locally 
fluctuating due to minor modifications of the experimental setup can be easily treated by a careful 
modelling of the variance/covariance matrix. 

We are studying extension of this method to Poisson uncertainties, and to the evaluation of the resolu- 
tion limits of a fluorescence TCSCP apparatus |7|]. The ultimate goal is to obtain consistent uncertainty 
estimation for lifetimes recovered from fluorescence spectra analysis. 

An alternative treatment is to model the IRF by a function, which parameters pdf 's are estimated by a 
bayesian analysis 

p(m) f 
p(m|s,Rs,h,R/j) = ^-— / dphP{s\m,Rs,Ph)piPh\K'Rh)- 
P(s) J 
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